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Abstract. Wc study the nonlinear dynamics of a warped or twisted accretion 
disc, in which the viscosity coefficients are assumed to be locally proportional 
to the rotational velocity prescription). Using asymptotic methods for thin 
discs, dynamical equations of the disc are obtained in a warped spherical polar 
coordinates. These equations are solved by the method of the separation of the 
variables. This analysis constitutes an analogous study of the nonlinear theory 
of an alpha model warped disc which has been studied by Ogilvie (1999). We 
have compared our results with Ogilvie's analysis. The dynamical behaviors 
of these models have also been discussed. Our results show that different vis- 
cosity prescriptions and magnitudes (a and /3 prescriptions) affect dynamics 
of a warped accretion disc. Therefore it can be important in determining the 
viscosity law even for a warped disc. 
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1. Introduction 

The development of studies about accretion discs during present century is an illustration 
of the grow of interest of researches. In more recent times a good deal of attention has 
been devoted to studies of the warped accretion discs. This kind of discs has been observed 
in a wide variety of astronomical objects from young stellar objects. X-ray binary stars 
to active galactic nuclei. For example, the stability of the super orbital periodicity in the 
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neutron star XRBs Cyg X-2, LMC X-4 and Her X-1 (Clarkson et al. 2003) is attributed 
to a warped disc. The outflows and illumination patterns from the central engine of a 
Seyfcrt Galaxy (Grccnhill et al. 2003) provide direct reasoning to warped accretion discs. 
Also, observations of FUSE (Far Ultraviolet Spectroscopic Explorer) and HST (Hubble 
Space Telescope) ultraviolet of the low-inclination, nova-like Cataclysmic variable RW Sex 
presented evidences of warped discs (Prinja et al. 2003). The precession of warped discs 
in magnetized stars (T Tauri stars, white dwarf or neutron stars) discussed by Pfeiffer & 
Lai (2004) that can be concluded from magnetic torques due to the interaction between 
the stellar field and the induced electric currents in the disc. A photometric study of the 
SW Sex star, PX And (Boffin et al. 2003) reveals a precessing disc possibly warped. Many 
different systems, including young stellar objects and X-ray binaries, display a different 
precession period or radiation flux variabilities that it may be concluded from a warped 
accretion disc. Also, Her X-1 (Tananbaum et al. 1972; Katz 1973; Roberts 1974; Still 
& Boyd 2004) and the hard X-ray component of the micro quasar GRS1915-I-105 (Rau, 
Greiner & McGoUough 2003) show another typical behavior of warped discs. 

During recent years, several driving mechanisms for disc warps have been suggested. 
A warp may be induced through instabilities that may be due to viscous torque (Pringle 
1992), bending waves and viscous torques (Ogilvie 1999), resonant tidal interactions 
(Lubow 1992; Lubow &: Ogilvie 2000); irradiation-driven wind torques (Schandl & Meyer 
1994); radiation torques (Pringle 1996); wind torques via Kelvin-Helmholts instability 
(Quillen 2001); magnetic torques (Lai 1999) and magnetic-induced electric current torques 
(Lai 2003). 

Most theoretical studies of the accretion discs are based on the concept of a real 
fluid and the equations of magnetohydrodynamics. Adjacent layers of a moving real fluid 
experience tangential forces (shearing stresses) as well as normal forces (pressures). It 
seems that the viscous forces are the main physical agent in any accretion disc theory. 
Therefore, the viscous forces are parameterized. A way of doing this is through the 
dimensionless parameter a introduced by Shakura & Sunyaev (1973). It describes angular 
momentum transport in accretion discs. The mechanisms of the angular momentum 
transport can be described by the magnetic or the hydrodynamic instabilities. There 
is another prescription through the dimensionless parameter /3 introduced by Lynden- 
Bell & Pringle (1974), in which the viscosity coefficient is proportional to the rotational 
velocity. Analyzing the turbulent flows between the coaxial cylinders, Richard & Zahn 
(1999) investigated /3- prescription for turbulent viscosity and found that the parameter f3 
is of the order of 10~^. Furthermore, in their model, the viscosity is very small comparing 
with the viscosity in a- prescription. Consequently, the forces due to viscous friction are 
very small compared with the remaining forces (gravity and pressure forces). So, the 
Reynolds numbers are very large, because of very low viscosity of the fluid. Accretion 
disc model with /?- prescription has been studied by Hure, Richard & Zahn (2001). They 
argued that the beta model for viscosity prescription is applicable for analyzing the steady 
state structure of the Keplerian accretion discs and it yields somewhat different results 
compared to the classical a-viscosity introduced by Shakura & Sunyaev (1973). Study 
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on hydrodynamic viscosity and selfgravitation in non-warped accretion discs (/? model) 
has been done by Duschl et al.(2000). They showed that /?-discs can explain the observed 
spectra of protoplanctary discs and yield a natural solution to an inconsistency in the 
a-disc models if the mass of the disc is large enough for self-gravity to play a role and 
in the limit of low mass, hydrodynamic turbulence will result by a model. Turbulence 
induced by the horizontal and vertical shear has been studied by Mathis ct al. (2004). 
They have presented a new prescription, the /3-viscosity, for the horizontal component 
of the turbulent diffusivity due to the differential rotation in latitude. They generalized 
(3 prescription (Richard & Zahn 1999) in the stationery limit, advcction and diffusion 
balance each other. Because that prescription (Richard & Zahn 1999) has been established 
in the case of maximum differential rotation and so its validity must be verified to milder 
shear rates. They have examined various prescriptions with their work. Their prescription 
yields a better agreement with the observations, but one can hardly consider it as the 
final answer, especially for extreme differential rotation. 

In general, what we can find from previous works on accretion discs is the effects 
of various viscosity prescriptions on the structure of discs and under conditions, their 
results arc compatible to some prescriptions. However, almost all previous studies of 
/3-discs dedicated to non-warped accretion discs. 

In order to investigate the non-linear dynamics of warped accretion discs and whether 
selecting of the model affects the forms of the equations governing a warped viscous disc, 
we applied Ogilvie's method (1999, hereafter OG) as considered in a a theory and check 
this new prescription for a thin viscous disc. Our study includes numerical solutions 
of the Navier-Stokes equations. First, the problem is reduced to a so-called singular 
perturbation which is then solved by the method of matched asymptotic expansions. 

The structure of the paper is as the following. In the next section, we explain gen- 
eral formulation. Analysis of the problem is presented in section 3. Then we show that 
the equations can be solved using the method of the separation of the variables and the 
numerical solutions are discussed in section 4. We compare the warped a and /3-discs in 
section 5 and then we conclude in section 6. 



2. General Formulation 

In order to construct a model for a warped accretion disc, we start by writing hydrody- 
namic equations. Then, we use the appropriate forms of these fundamental equations in 
warped spherical polar coordinates {r,9,(f)) (OG). Although many physical agents such 
as magnetic fields and radiative processes are playing significant roles in the dynamics of 
the discs, we neglect all those complex phenomena in order to understand the dynamics 
of an accretion disc with /3-prescription for viscosity via semi-analytical methods. Also 
the self-gravity of the disc and interaction of the stellar magnetic field with the disc are 
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neglected. The fundamental governing equations are the continuity, 

Dp = -pV ■ u, (1) 

the adiabatic condition, 

Dp = -TpV ■ u, (2) 

and finally the equation of motion, 

pDu = -Vp - pV$ + V • [/iVu + /i(Vu)'^] + V [(^ib - |/i)V • u] , (3) 

where u, p, p and $ are the absolute velocity, the density, the pressure and the external 
gravitational potential, respectively. Also, T is the adiabatic exponent and the shear and 
the bulk viscosities are denoted by /i and /ib- The symbol D denotes the Lagrangian time 
derivative operator, 

^ = {dt)r,e,<i> + vrdr + -de + ^^d^. (4) 
r rsmd 

Note that the components (wr, ug,u^) are the absolute velocity components, i.e. u is the 
velocity as measured in the incrtial frame. But the additional motion with respect to the 
warped coordinate system is described by the relative velocity v. The exact relationship 
between these two velocities has been found by OG. 

Warped spherical polar coordinates (r, 0, </>) is defined so that on each sphere r = 
constant, one can define the usual angular coordinates {9, (j)), but with respect to an axis 
that is tilted to a point in the direction of unit vector £{r,t). This tilt vector can be 
described using the Euler angles j3E{r,t) and 7(r, t): 

£ = sin /3b cos 7 ex + sin /3b sin 7 + cos (SeGz- (5) 

For the viscosity coefficients, we arc using the /?- prescription (Lyndcn-BcU and Pringlc 
1974) rather than the usual a- prescription which has been used by OG. Thus, the 
viscosity coefficients are assumed to be locally proportional to the rotational velocity, 

M = Pr^^P, (6) 
Atb = PhT^^P, (7) 

where the dimensionless coefficients /? and /3b can be considered as functions of the radius. 
We will show that (3- prescription for the viscosity coefficients lead to significant changes in 
the equations describing a warped disc. Also the fiuid is assumed to be locally polytropic, 

p = kp^. (8) 

where r(r) is a prescribed function of radius. 



Now we can consider a thin disc in a spherically potential <i?(r), in which the small 
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parameter e is a characteristic value of the local angular semi-thickness of the disc. Using 
this small parameter, it is possible to study the structure of a warped disc by asymptotic 
expansion method. In thin disc approximation, one can assume that the disc matter 
lies close to = 7r/2. To resolve the internal structure of the disc, introduce the scaled 
dimensionless vertical coordinate Q, 

^ = f-<, (9) 

and the slow time coordinate, 

T = eh. (10) 
For the density and pressure, introduce the scalings 

p{r,6,<l>,t) = e'[p<y{r,4>,C,T) + epr{r,4>,(:,T) + 0{e^)\, (11) 
p{r,e,cl,,t) = e«+2[po(r,0,C,T) + epi(r,0,C,T) + O(e2)], (12) 

and for the relative velocities, 

Vr{r,e,<(>,t) = et;,i(r,0,C,r) + eV2(r,0,C,r) + O(e3), (13) 
ve{T,e,<t>,t) = evei{r,(i),Q,T)+e\e2{T,(i),Q,T) + 0{e^), (14) 
v^{r,e,(t),t) = rf2(r) sin 61 + ev0i(r,0,C,T) + e%2(r,0,C,r) + O(e^). (15) 

Finally, for the viscosities, assume 

Kr,0,(/),t) = e'+^[MrA,C,T) + etMr,^,C,T)+Oie^)], (16) 
tib{r,e,<jj,t) = e«+2 [^bo(r,(^,C,I^) +eA*bi(^, r) + 0(e2)] . (17) 

where .s is a parameter which should be positive if the self-gravitation of the disc is to 
be negligible. The equations of fluid dynamics were derived in warped spherical polar 
coordinates by OG. He reduced them by means of above asymptotic expansions for a 
thin disc and then divided them into two sets. Set A, which determines the intermediate 
velocities, consists of five coupled non-linear partial differential equations (PDEs) in two 
dimensions ((6, C) and seven dependent variables {pQ,p(). I-J.q. i-j.ba,Vri,iJoi,v^i}. Also Set 
B, which determines the slow velocities, contains a set of five linear PDEs for the higher- 
order quantities {pi,pi, fii, p.bi,Vr2,V02,Vii,2}, including coefficients that depend upon the 
solutions of Set A and their radial derivatives. In the present work, we adopted Set A 
and B for our analysis (see OG, for details of deriving the expansions of Set A and Set B). 



3. Analysis 



It is very unlikely that the equations of Set A can be solved analytically. Numerical 
approach is a convenient way. Fortunately, we can transform the partial differential 
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equations into a set of ordinary differential equations (ODEs) using the method of sep- 
aration of variables. To achieve this, at first we introduce the following forms for the 
variables: 

= r^^^[h{<ly~x)~\f2{<ly~x)e], (18) 

Vri = rl]/3(</)-x)C, (19) 
vei = rl][/4(</)-x)C + .94(<^-x)C'] , (20) 
v^i+rVrW cofi 13 E = rf] [/5((/)-x)C + 55('/>-x)C^] • (21) 

where /iq = p3Y^ is the enthalpy. Using the relation (|18|l . one can drive the upper 
surface of the disc as 

e = 2/i(0-x)/2"'(0-x) (22) 

Introducing by the dimensionless functions /i, . . . , /s and g^, g5, one can obtain the cou- 
pled sets of the non- linear ODEs of first order for the given equations 





= (r-i)/4(0)/i(0), 




(23) 




= (r + i),ui<p)f2W - 6(r - i)g,{^),h{^), 




(24) 


I'M 


= fMfm + 2/5(0) + [/2(</>) - 6(/3b + \P)gm 


] V COS0 


(25) 


m) 


= -fM\^\ cos + 2/3(<^)|V'| sin <^ + /4(0) [fi{4>) - 


|-/3(0)|V'|COS 0] 


+ 1-/2(0) 




+6(/3b + i/3).g4(0) + 6/334(0) (1 + cos" 4>), 




(26) 




= .94(0)/3(0)|V'| COS0 + Agi{(t>)fi{(t)), 




(27) 


m) 


= UWfsW - i^iV3(0) +6/3g5(0)(l + |Vpcos2<^ 


'), 


(28) 




- 3/4(0)55(0) +54(0)/5(0), 




(29) 


where these functions are subject to periodic boundary conditions /„(27r) 


= /„(0) and 



f/„(27r) = 5n(0). The epicyclic frequency K(r) is defined by 



K 



•2 _ /io2 



+ 2rnn', (so) 



and the dimensionless epicyclic frequency is k = n/n. Meanwhile the amplitude of the 
warp is defined as 

IV'I = r\^^l (31) 

hence, using the tilt vector, we have 

1/, = |V;|e'>^ =r(/3^+i7'sin/3E), (32) 

that it is a dimensionless complex variable.^ Moreover, by the definition of /g (see. 
Appendix B), one can drive 

/a0) = -2/4(0)/6(0) + 9(r - l)/l(0)/2(0)-'g4(0)/4(0). (33) 



^Throughout this paper, a prime and a dot for Q{r),f3E and 7 imply radial and time derivatives, 
respectively. 
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Also in addition to these, we find the foUowing combinations 



/2 (0)54(0) 
54(0)/3(0) 



if r 7^ 1/3, 
-255(0) all r. 



(34) 
(35) 



All information about a warped disc with (3- prescription can be obtained by solving 
equations H23l) - (|33|1 numerically. Before presenting the results of numerical integration, 
we discuss about evolutionary equations of a warped /3- disc and its relations with solu- 
tions described by functions /„ and g^- In particular, it is very important to see whether 
the equations of Pringle (1992) can be derived from the three dimensional fluid equa- 
tions. Pringle (1992) developed an approach for deriving the equations of a warped disc, 
without reference to the detailed internal fluid equations. OG showed the impossibility 
of deriving the angular momentum equation of Pringle (1992) from the basic equations of 
fluid dynamics. In fact, one should allow for a more general form of the torque between 
neighboring rings. Thus, OG introduced three dimensionless coefficients Qi, Q2 and Qa, 
which depend on physical quantities of the disc. These coefficients enable us to discuss 
about relative importance of torques between neighboring rings in a warped a- disc. We 
are following similar approach for describing a warped (3- disc. However, as we will see, in 
a warped [3- disc there are various torques and so we should introduce extra dimensionless 
coefficients. Using the equations that can be extracted from Set B in OG and substituting 
the defined quantities {/lo: 'f^ri, I'ei, I'^i}, we can propose the evolutionary equations for 
the warped (3- disc as (for details, see Appendix B), 



ot ' r or ^ ' 



ld_ 

r dr 

+ - — 

r dr 



{Qiir^n^e) + 
{Q[j:r'^n^e) 



ld_ 

r dr 

+ - — 

r dr 



dr 



Q'^T^r^n^ 



m 

dr 



r or or 
r dr \ dr 



(36) 



for the angular momentum, and 



d 



1 d 



dr r dr 



{Qilr^n^) - Q2lr^^'^ 



de 

dr 



1 d 



r dr 



m 

dr 



for the component of angular momentum parallel to £, and 



(37) 



or r or \ or 



dl 
dr 



r or \ or 



or r or \ or 



d£ 
dr 



r dr \ dr 



(38) 
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for the tilt vector. 

Thus, we can obtain coefficients Qn and QJ^ {i.e.,n = 1,4} as (Appendix B), 

Qi = (/6(-/3/5 + 3/355|V'|cOS</))), (39) 

Q[ = {~^{4~k^)P + pU\ij\cos^), (40) 
Q4 = ^(e'^/6[/3-i/3(/4 + /3|^|cos0)+3i/354|V'|cos0]), (41) 

Q'i = ^(e'^[i/?(/4 + /3|V'|cos0)|V'|cos0-i/?/3-i/3|V|sin0]), (42) 
where Q4 and Q'^ defined as Q2 + iQa and Q2 + iQ'^. 



The dimensionless coefficients Q„ and Q'^ can be evaluated from the solutions of Set 
A (see OG). We find that these coefficients depend on selecting the model, the amplitude 
of the warp, the rotation law and the shear viscosity. For a warped /?-disc as we consider 
here, we find three kinds of the internal viscous torques, (see equation l|3()|ll: 

G{r,t) = {QiI+Q[j:r^yn^l+iQ2X+Q'2^r^yn^^ 

or 

8f 
or 

(43) 

From a mathematical point of view, the equation H3fci|) constitutes the prototype for 
a parabolic partial differential equation and can be thought as an advection-diffusion- 
dispersion equation in the non-linear regime. So for a flat disc (^ = 0), the evolution 
equation reduces to a standard disc diffusion equation. 

According to the above equation for G(r,t), the first term on the right hand side gives 
a contribution to G which is in the local direction £. So, the effective advection coefficients, 
Qi and Q'l, represent viscous torques on each ring in the disc due to differential rotation 
within the disc plane. Thus, the rings tend to rotate in direction £. 

Considering neighboring rings in directions £ and £ + A£, the second term on the right 
hand side gives a contribution to G which is in the |^ direction. Therefore, the effective 
diffusion coefficients, Q2 and Q2, represent viscous torques due to the interaction between 
neighboring rings in the disc, in order to flatten the disc. 

The last term on the right-hand side is a dispersion one. So the effective dispersion 
coefficients Q3 and Q3 demonstrate torques which are perpendicular to £ and ^ , respec- 
tively. In this case, each ring in the disc experiences torques tending to make the ring 
precess if it is not aligned with its neighbors. Then we would expect to generate wave 
motions in the disc so that the warp propagates as a dispersive wave. 
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4. Numerical solutions 



One can expand the dimensionless functions /„ {n = 1,...,6}, and 55 as power 
series of the amplitude of the warp to determine the dynamics to any desired order. 
In Appendix A, truncated Taylor series are presented including the coefficients Qn and 
Q^. These scries enable us to start numerical integration. From the expansions of the 
dimensionless functions /„ and Qn and their relations to coefHcients Qn and Q'^, we found 
truncated Taylor series only for coefficients Q[,Q2 and Q3. Meanwhile, as we discuss in 
Appendix A, our solutions are restricted to the nearly Kepler ian-discs. So we consider in 
our calculations for = 0.99 for all values of (3. 

In order to be able to compare our results with the warped a- disc, we adopt coef- 
ficients Qi and Qa = Q2 + iQs from OG. To distinct between these two coefficients, we 
label them by superscript a. Now the set of ODEs in the previous section can be solved 
numerically and we may obtain the coefficients Q'^ and Q„. In addition, we should note 
that solutions must satisfy periodic conditions for functions /„ and gn as well as (/g) = 1 
(see, Appendix B). Thus, the derived solutions make an interesting dynamical description 
of a warped /3- disc. 



4.1 The warped /3-disc 



First, we investigate a non-Keplerian disc without viscosity. We consider F = 5/3 and 
(3 = (3b = 0. The dimensionless functions /i, /2, /s, /e and 54 are even, but /s, fi and 55 
are odd. So the only non vanishing coeSicient is Q3. So, 



Qs = 77T(/6[/3Sin0-/3COS(/)(/4 + /3|V'|cOS0)+3/3c/4|'i/'|cOS^0]) 1/3=0, 



(/e [/s sin ^ - /a cos (^(/4 + /a | V' I cos 0)] ) . 



(44) 



A contour plot of coefficient Q3, is shown in Figure 1. 




Figure l.The square of dimensionless epicyclic frequency versus the amplitude of the warp 
for Contour plots of the coefficients Q3 (solid line) and (dash line) for an inviscid disc 
with r — 5/3. The condition = 1 separates solutions into two parts. For < 1, we 
have Qs, > 0; and for k/^ > 1, we see Qs, Qg < 0. 

For the case of nearly Keplerian disc with viscosity, we consider — 0.99, F = 5/3 
and Pb = 0. Figure 2, shows /3 vs. \ip\ and contours of all coefhcients. Note that all 
the plots in Figure 2; parts (a), (c) and (e) show coefficients and parts (b), (d) and 
(f) show coefficients Qn ■ The solutions can be calculated for large values of 1-01. For 
evaluating the coefficients Qi and Q[, we used equations H39() and H40() . Figures 2a and 
2b show contour of plots of the coefficients Q[ and Qi respectively. For reasonably small 
values of there is a good agreement with the truncated Taylor series for Q[. 

Figures 2c and 2d show contour plots of the coefficients Q2 and Q2 respectively. We 
used equations H41|) and H42|l for evaluating the coefficients Q2 and Q'2, respectively. For 
reasonably small values of jV'l there is a good agreement with the truncated Taylor series 
for Q'2. 



Figures 2e and 2f show contours of the coefficients Qg and Q3 respectively. We used 
equations H41|l and (|42l) for evaluating the coefficients Q3 and Qg, respectively. Also Q3 
is much smaller in magnitude than Q3. For reasonably small values of I"!/)! there is a good 
agreement with the truncated Taylor series for Q3. 
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Figure2. The dimensionless viscosity parameter versus the amplitude of the warp for 
Contour plots of the coefficients; (a)Q'i, (b)Qi, (c)Q'2, (d)Q2, (e)Q'^ and (fJQs for a 
viscous, nearly Keplerian disc (fi^ — 0.99)with T ~ 5/3 and (3h — 0. 
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4.2 The warped a-disc 

First, as an illustrative case, we investigate a non-Keplerian disc without viscosity. We 

consider F = 5/3 and a = = 0. According to obtained functions fi, f2, fs, fe by OG, 
the only non vanishing coefficient is Q^. So, 

Qs = ^(/6[/3Sin(/.-/3COS(/)(/4 + /3|V|cOS<^)]). (45) 

If we plot n'^/Cl'^ versus we get the contours indicated in Figure 1 for Q^. We selected 
a given interval for |^|. 

For the case of nearly Keplerian disc with viscosity, we consider = 0.99, F = 5/3 
and ab — 0. Figure 3 shows a vs. \iJj\ and contours of all coefficients. Note that all the 
plots in Figure 3; parts (a), (b) and (c) show coefficients Q^. We adopt this coefficients 
from OG, 

Qi = (/6[^(4-«')a/2-/3/5+a/2/5|V'|cos<^]>. (46) 

•32 = ^</6[(/4 + /3|V'|cos</.)(l + /3|^|sin</.)+a/2/3|V|sin</. 

-a/2(/4 + /sIV'l cos(^)|V|' cos(/.sin<^ + a/ajVP sin^ ^]). (47) 

Qs = ■j;^(/6[/3Sin^-/3COS^!>(/4 + /3|V'|cOS(?i) + a/2(/4 + /3|V'|cOS(^)|V'|cOS^^!> 

-othhcoscj) - a/2|-0| sin./) cos (/>]). (48) 

5. Comparing the warped a and /3-discs 

5.1 Inviscid, non-Keplerian disc 

Figure 1 shows a good agreement between two contour plots of (^3 and Qg. Also they 

show good agreements with tnmcated Taylor series for small values oi As Figure 
1 demonstrates the condition = 1 separates solutions into two parts. For < 1, 
we have Q3,Q^ > 0; and for k"^ > 1, we see Q3,Q^ < 0. In the case < 1 and for 
small values of \tp\, the physical solutions of /? and a models exist only for \tp\ < 0.35 and 
[■^l < 0.4 respectively. The solutions terminates when 54 = 0, /2 ^ 1 for (3 model because 
we couldn't obtain the numerical solution. Whereas the solutions terminate for a model 
when /2 = 0. Hence, it seems that the stability condition establishes for such discs for all 
values > in small values of In the other word, under the condition > 0, since 
the epicyclic frequency has the relation with particle orbits, the displacements oscillate 
about a fixed mean position and the circular orbit is stable to small perturbations. 
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Figure 3. The dimensionless viscosity parameter versus the amplitude of the warp for Con- 
tour plots of the coefficients; (a)Qi, (b)Q2 and (c)Q'^ for a viscous, nearly Keplerian disc 
(k^ = 0.99 )with r = 5/3 and = 0. 
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5.2 Viscous, nearly Keplerian disc 

Now, we compare a and [3 models for a viscous, nearly Keplerian disc. If we compare the 
equation H3()|l with Pringle's equation H13U|I . the coefficients Q'^ and Qi represent viscous 
torques according to horizontal shear. Figures 2a and 2b show similar behaviors, the 
magnitudes of Q'l and Qi increase with increasing \'ip\ and this is similar to (Figure 
3a). Nevertheless, we see different treatments between two models when torques due to 
Q'l and Qi are considered. As mentioned before, they indicate viscous torques parallel to 
I and Q'l , Qi identified as the advection coefficients. Then comparing to a model, the 
advection of the warp occurs more significantly. On the other hand, the viscous stability 
condition implies that Qi, Qi < 0, so /3-discs comparing to a-discs tend to become stable 
more quickly. This is in agreement to Hure et al. (2001) based on viscous stability for (3- 
discs which is one of the important property of them. Thus, they do not tend to fragment 
in contrary to a-disc. 

Also, the coefficients Q'2 and Q2 represent viscous torques due to the vertical shear. 
Figures 2c and 3b show different behaviors. We see that with increasing j?/;!, increases 
whereas Q2 decreases. So, the results due to the torques of and Q2 don't tend to OG's 
results. As if, this model affects the vertical structure of a-disc so that viscous torques 
play more role to flatten /3-disc. Also comparing to the OG's analysis, the coefficients of 
Qg and Qg represent torques which lead to the dispersive wave-like propagation of the 
warp. Figures 2e, 2f and 3c show that Q3, and respectively, they treat in a similar 
manner so that with increasing they decrease. Although their typical treatments are 
similar, is typically much larger than Q'^ and Q'^ coefficients especially in small \'ip\. 
So we expect to generate wave motion more in small \'ip\ for /?-disc so that it loses its 
importance for large \'ip\. 

6. Summary and Discussion 

We have presented an analysis of the non-linear dynamics of a warped accretion disc using 
/3-prescription. Using basic equations of the fluid dynamics in a warped spherical polar 
coordinates, we obtained the general equations that describe a warped disc. We have also 
employed the method of matched asymptotic expansions for thin discs to derive a set of 
coupled PDEs which govern the dynamics of the system. Then we solved the equations by 
the method of the separation of variables in order to extracting the equations governing 
the warp in their simplest forms. The non-linear dynamics of a warped /?- disc under a 
differentially rotation field in the presence of a spherically symmetric external potential 
presents a more complicated behavior comparing to the case in which a warped a- disc is 
considered (OG). Moreover, it is more complicated than the case of the linear theory of 
Papaloizou & Pringle (1983) to allow an arbitrary rotation law. In our study, there are 
six coefficients. They have been determined numerically and analytically, by solving a 
set of ODEs. Using truncated Taylor series in the amplitude of the warp, the coefficients 
Q'l J Q'2 ^-nd Q3 are calculated analytically. Figures 1-3, show their contour plots for two 
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cases inviscid-Keplerian and viscous-near ly Keplerian disc when F = 5/3 and /Sh = 0. 
We found that these coefficients depend on selecting of the model and the shear viscos- 
ity. Considering the equations (jSHJ ■ ■ • ^ one can easily note that the coefficients also 
depend on F, and the value of Our results show that the equations governing a 
warp viscous disc depend on the parameters of the model. The results can be compared 
with the studies of Pringle (1992) and OG as follows: 

1) A comparison with the notation of Pringle (1992), there are four different viscous 
torques due to the interaction between neighboring rings in the disc (see, Appendix B). 
Therefore, these torques may be explained by four coefficients. Comparing to OG, we 
will find two extra coefficients Q'l Q'2 ^^r the viscous torques. To understand the 
meaning of these coefficients, we may deduce following equations by comparing equations 

(innji and (nnnii, 

QiXr! + Q'lSr^r! ^ I'lS^T^' (49) 

amr 

Q2in + Q'2j:r^n < — > ^^2^. (50) 

We see that the numerical evaluation of the coefficients Qi,Q'i and Q2, Q'2 represents the 
qualitative behavior of vi and 1/2- We may conclude that i^i and 1^2 depend on the value 
of ^. 

2) The evolutionary equations for the warped disc show that this scheme is a general- 
ization of the form suggested by OG. A warped accretion disc can be studied in detail 
by the coefficients Q3 and in (3 model and only the coefficient Q" in a model (OG). 
These coefficients demonstrate torques tending to make the ring precess if it is misaligned 
with its neighbors. Comparing Figures 2e and 3c show that the direction of the torques 
is negative for both coefficients Q'^, Q". These torques lead to the dispersive wave- like 
propagation of the warp. 

3) Our results show that different viscosity prescriptions and magnitudes (a and (3 pre- 
scriptions) affect dynamics of a warped accretion disc. Therefore, it can be important in 
determining the viscosity law even for a warped disc. 

Some further work still remains to study the dynamics of a viscous warped accretion disc 
using /3-prescription: 

a) The dynamics of the warped accretion disc may be studied by considering the mag- 
netic effects. Then we have to modify the equations by including magnetic field of the 
disc. Thus, one should introduce a suitable model for the geometry of the magnetic field 
components. Also, one can extend our analysis to the case that self-gravity of the warped 
discs is important. 

b) In this paper, we neglected the thermal and the radiative effects, however, these phys- 
ical processes are playing important roles in the dynamics of the discs. Therefore, the 



16 



J. Ghanbari , M. Shadmehri and F. Salehi 



behavior of the dynamics of these discs needs to be studied further. 

c)As we discussed, the expansions fail only when = 1 for all values of f3. This is 
a resonant case, which can't be described using this method. Therefore, it is interesting 
to study the non-linear dynamics of the resonant case for a viscous Keplerian (or nearly 
Keplerian) disc. 
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A. Truncated non-linear equations 

In this Appendix we find truncated Taylor series of the dimensionless functions Jn{n — 
1, . . . , 6}, 54, 55 and the coefficients Qi, Q[, Q4 and Q'^ in terms of powers of {ipl. Thus, 



/l(0) 


= /lO + IV^lVl2(0)+O(|V^|4), 




(51) 


/2(0) 


= /2O + |^|V22(0)+O(|^|4), 




(52) 




= ^|/3l(^) + |V;|2/32(</') + |V'lV33(0)- 




(53) 




= VP/42(0)+O(|^|"), 




(54) 


54(0) 


= ^1341(0) + |Vl'343(0)+O(|^|'), 




(55) 




= ^1/51(0) + |V'P/52(0) + |V'l'/53(0)- 




(56) 


55(0) 


= ^Pg52(0) + |V^Pg53(0) + O(|vn, 




(57) 


M<t>) 


= /60 + |^P/62(^) +0(1^-1')- 




(58) 



Note all terms that are scaled with \tp\'^ indicate an unwarped disc. 
A.l Zeroth-order solution 

For an unwarped disc, equation (|26|l at 0{\ip\'^) gives 

/20 = 1. 

Also equation at 0(|V'|*') yields 

and since (/g) = 1 by definition (see. Appendix B), then 

/60 = 1. 

A. 2 First-order solution 

The vertical velocity at first order is obtained using equation (|27|l at 0{\ip\), 

9'aM - 0, (62) 

Hence 

341 = Cei- (63) 



(59) 
(60) 
(61) 



18 



J. Ghanbari , M. Shadmehri and F. Salehi 



Also the horizontal velocities at first order are obtained by equations H25|l and (|28|l at 

o(IVI), 



f'^M-V^oi{4>) = CO 



Therefore, 



fsiW = Cri cos (/i + S'ri sin 0, 
/5i(</>) = C^i cost/) + 5^1 sin (/). 

We can express the solutions using a complex notation as Z = C + iS", 



" -i -2' 




Zrl 




■f" 














(64) 
(65) 



(66) 
(67) 



(68) 



The determinant of this matrix is —(1 — k^), and therefore, there isn't any solution for 
the Keplerian-disc (k^ = 1). Hence,we have 



Zd,\ — 



2(1 -k2)- 



(69) 
(70) 



A. 3 Second-order solution 



The horizontal velocities at second order are obtained by equations (|25|l . 1)28(1 and H29|l 

at 0{W) 



/32(</')- 2/52(0) 



-6(/3b + \li)g^\ COS0, 



/52('/') + 5^'/32(0) = 6/3.g52(0), 
552('/') = 341/51 ((/>)• 

Equation yields a solution as 

552((/>) = 1^01 sin./) 
while equations H71|l and (|72|l have the solutions as follows 

/32(</>) = Cr2 COS </i + 5*^2 sin 

/52(</>) = 6*02 COS (/i + 5^2 sin (/), 



(71) 
(72) 
(73) 



(74) 



(75) 
(76) 
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with 



so the solution is 



" -i -2' 




Zrl 













Zr2 — SiCei 



-6(/3b + \P)Cei 



^02 



;i2(/3b + i/3) + 2/3Z^i 



^2 - 1 



(77) 

(78) 
(79) 



Also the enthalpy and vertical velocity at the second order are obtained by equation 

at 0(|V^|2), 



f[2W = (r-l)/42(<^)/l0, 

equation at O(|-0P), 

/22('^) = (r + l)/42(0)-6(r-l)/io54i, 
and equation at 0{\ip\'^) 

/42(0) - -/^i(0) cos 0+ 2/31 (0) sin 0-/22(0). 
Combining the last two equations, gives 

/42(</') + (r + l)/42(0) = 35,isin20. 
It has a solution as the following form 

/42(0) Ce2 cos 20 + 5^2 sin 20, 

with 



(3 - r) 




Ce2 







-(3-r) 




. ^S2 _ 




3S^\ 



3i 

r - 3 



Sri- 



In a complex notation, the solution is 

.^612 

It then follows that 

/l2(0) = -i(r-l)/io5e2COs20, 

/22(0) = -i(r+l)5e2cos20+-^5,i, 



r-3 



/62(0) = 5*92 cos 20. 



(80) 

(81) 

(82) 

(83) 

(84) 
(85) 

(86) 

(87) 
(88) 
(89) 



20 



J. Ghanbari , M. Shadmehri and F. Salehi 



A. 4 Third-order solution 

The horizontal velocities at third order are obtained by equations (j25(l . H28|l and H29|l at 



fsM-V^sW = /42(</')/3l(0) + /22(</')cOS0, 
+ 5^'/33(0) = /42(0)/5l(</')+6/3<?53('/'), 
553 ('^) = 341/52(0)- 

Equation H92(l yields a solution as 

953{(l>) = 6*91^02 sine/) 

while equations H9U|) and (|91|l have the solutions as follows 

733(0) = C'rs cos (j) + Sr3 sin + {m = 3 terms}, 
/53(0) = COS + sin {m = Sterms}. 



hence, one can obtain 



k2 - 1 
1 



- 1 



r r + 1 1 ~ ' 

— — ^S'ri H — Se2 ~ ■^Sg2Sri + Se2Z^i + 12(iCeiZ^2 

r + 1 1 

Srl Se2 + 



25*92^01 + 6/3C*eiZ02 — 2" V T - 3 



(90) 
(91) 
(92) 



(93) 



(94) 
(95) 



(96) 
97) 



The vertical velocity at third order is obtained by equation 127(1 at 0(|V'P) 

543(0) = 4(?4i/42(0)+ 541/31(0)543(0) cos (98) 
it yields a solution as 



543(0) — —2S02C01 — -SriCei- 



(99) 



A. 5 Evaluation of the coefficients 

With respect to the functions /i-.-/6,54 and 35, we find truncated Taylor series for the 
coefficients Q[ and as: 

Q'l = Q'io + l'/'l'Qi2 + o(|^|^), (100) 

Q'a = O4O+ 1^1^41 + IV'l'Q42+0(|^n. (101) 
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At zeroth order, we obtain 



and 



0'io=^(« -4)/3, (102) 
g^o = -i/3(e'^(/3iW+sin<^)) 



At first order, 



Q^i = -i/3(e'^/32(</.)) 



2 

= -3/3C.,^^ii±^^p-^^^^^. (104) 



At second order. 



(/3b + i/3)(K2-l) + /3K2 
1 



Q'l2 — 2^^"^^ 

/3 K 



^•2 



4 (1 - k2) 



(105) 



and 



Q42 = i/3(e'^/42(<^)cos0 + /3i(</.)cos2(/,-/33((/.))> 



^ -{4r(K2 - 1) - 3(r + 1)(k2 - 1) - 6 - 6k2(k2 - 1) 



8(k2 - l)3(r-3) 

+144(r - 2,)pCl^[R\f^ - l)(/?b + i/3) - PA). (106) 

From the coefficients Qi and (^4, truncated Taylor series were found only for the 
latter, 

Qi = g4o + |V'|Q4i + |V'l'Q42 + 0(1^-1^). (107) 

At zeroth order, 

(340 = (e'*/3i (</>)) 
= 2^'-i 



2(1 -k2) 
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At first order, 



Q41 = (e'^(/32W+3i/354icos<^)) 

3i^ 3(/3b + l/3)(;^^-l)^ + /?K^ + /3(;^^-l)3 



At second order, 



Q42 = (e'^/33(<A) - i/3l(</')/42(<^) - cos<^)) 

= -iPQ42- (110) 



B. Evaluation of the general form of the angular momentum 

equation 

We wish to drive the general form of the angular momentum equation. At first we adopt 
the equations from OG. ^ 



or 



J {'^^'^ [/"o^^fi' - Por^Vri{v4,\ + rVrW COS (3e) 
+/zor^(/3^ cos(^ + j'sitiPe sin 0)9^(^01 + rVriY cos f3E)] )rdC 
-J {Po b'01 + rvri {P'e cos (f) + -f' sin Pe sin (p)] [rQ( + rvri {P'e sin </> - 7' sin Pe cos 
—PqWe sin 4> ^ 1' sin (3e cos <^)d(^Vri 

+/xor(/3^ cos 4> + y sihPe sin 4>){I3'e sin0 — 7' sin (3e cos 4>)dc^ [vqi 
+rvri {P'e cos + sin Pe sin ^)] — /xor^n(/3^ sin </> — 7' sin /3b cos (/>)^ 



)rdC. 

(Ill) 



^They can be extracted from the Set B by integration. Note that, the operation (.) stands for 
azimuthally averaged quantities. The range of the integrations with respect to <t> and C, are from to 27r 

and —oo to oo, respectively. 
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and 

X [vei + rvri {P'e cos + 7' sin /3b sin 4>)] 

+i/ior"^ cos + 7' sin /?E sin 0)9^ [t;^! + ru,.i (/?^ cos + 7' sin sin 0)] 
-i^o7-^9<;Wri - iAfo?''*^^ sin - 7' sin /3s cos 0) } ) r dC + j (e'"^ { -po (r^f7)'Cwri 
-Po(w0i + rWri7'cos/3E) + to^i cos (/) + 7' sin/3^; sin0)] 

+ipO?'f^C(''^01 + ?'Wrl7'cOS;9£;) 

-\-\porVri{v^i + rvrij' cos /3b) (/3^ sine/) - 7' sin/3^; cos</)) 

+ rUri7'cos/3B) - i^o^^C/^B sin0 - 7' sin /3b cos(/)) 
X [rD,' + (/3b COS0 + 7' sin /3b sin(/))5^(u0i + ru,.i7' cos /3b)]} )rdC- 

(112) 

where 

Swr = J {poVr2 + PiVri)rd(, = J po r (113) 

here Wr(f', t) and I](r, t) are the mean radial velocity and the surface density, respectively. 
To proceed, we define 

/6(</>-x)-i/J (114) 
where T{r, t) is the second vertical moment of the density being defined through 

I = y PoC'^'dC (115) 

so that I is its azimuthal average, i.e. I =< I > and also the definition of /g requires 
< /g >= 1. Therefore, by Substituting the relations and three 

defined relations Q and (|114|l into the equations (|111|1 and pi2l) . the first integral 
in pil|) becomes 

+3{r^n^pg5fe\i^\cos^)li 

(116) 

and the second integral in (|lll|l becomes 

h = -(172/4/6)2:- (02/3/4/61^1 sin <^)X-(r2f]2;j/3|V'|sin0)E 
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+{r^VL^I3f4\ii\'^ cos (j> sin 4>)T, + [3gih\i>\'^ cos sin 0) J 
+(r202/3/3|i/>|3 cos^ 0sin<;i)E - {r^ ^'^ I3\^P\'^ sin^ 

- ( O Vs/e I I cos 0)X - ( r! VI /6 1 V' I ' sin cpcos^)l. (117) 

similarly, the first integral in H112(l [Note that, for any function F, we can write, 
X)) = {e'^'t'+^^F{4,)) = ^(e''^F(0))] becomes 

- i(a,+i7'cos/3B)[(e'Vl72(l-i/4-i/3|7/;|cos0)/3/6)l^ 
+i(e'Vo2/?/4|V;| cos 0)1]^ + 3i(e'^2f^2^g4/6|^| cos 
+i(e'Vf72/3/3|V.pcos2</>)S^-i(e'Vr!2/3(/3 + |^|sin^))l]^], (118) 
and finally for the second integral in (|112|l .we have 

I'2 = -(e'^y f^Va/e)!^ - {e'*^\hh + /sM^I cos^)/^)!^ + i(e'*f]2^5/6)l^ 

+i(e'*l]V3/5/6|V'| sin0)l^ - {e'\^V?(3h)^^^ - ^{e'^n' fihg^)!^^ 

-i/ei-^^f^^r^rj^/?!^! sin(^)S^ - i/e'^r^riVs/JlV-P sin0cos0)l]^ 
2 'IV'I IV'I 

J, 

-3i(e"*rj2g5/6/3|V'l^sin(/)COS0)X|;^. (119) 
in which /„ stands for /„(</>). 

Now, we attempt to arrange them in terms of coefficients Qi,Q2, Q'l, Q'2 a-nd (54 = 
Q2 + iQs, Q4 = Q2 + iQs- Then, we have 

h = iaJQir^f^^J+Qir^'f^^S], (120) 
r 

h = -Q2f^'2:|V'P-Q2^'f^'s|i^|2. (121) 



and also 



I[ = -{dr + W cos l3E){QAr^^'^Ii' + Q'ir'^^'^^i'), (122) 
r 

I2 = Qim^tp + Q[T.r^n^i:. (123) 



Therefore, it is possible to get the coefficients Qn and Q'^li.e., n = 1, . . . , 4} in terms 
of / and g. In the next step, we may write the following combinations, 

or 
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r or 



de 

dr 



de 

dr 



ot or 



T' + J' 



„de 19/ -3 -^dt \ on 

or r or \ or 

+ ^-^(Q,ir^n'ex^ 

r or \ or 



(124) 



de 

dr 



de ^ ld_ 
dr r dr \ ^ 
dt 
dr , 



1 d 



de 

dr 



--^ ( Q'^r^n^ex^ 

r dr ' 



(125) 



where we used relations ()ril|l. (|32|l . Also it has been considered that the disc matter 
will lie close to = n/2. So we could write £ in terms of spherical polar unit vectors 
(er,ee,e0), 



-eel 



(126) 



hence 



de 

— = — er(/3£; C0S(7!) + 7sin/?£; sini/)) + e0[/3£; sin(/) — 7sin/3£; sin0 cos(/)], (127) 



de 

dr 



-BriP^ cos(p + 7' sin/3^; sini^) + e0[/?^ sin0 — 7' sin (3e sin0cos(/)]. (128) 



To compare the present work with Pringle's work (1992), we present the following 
equations (Pringle,1992), 



for the surface density E(r, t), and 

ot r or r or \ dr 



1 a A ^ 3„ar 



(129) 



(130) 
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for the angular momentum L = YiR^Vll, in the absence of external torques. Here vi 
and V2 are the viscosity corresponding to the azimuthal and vertical shears, respectively. 
From these equations, one can derive 



or r Or 



dr 



dr 



(131) 



for the component of angular momentum parallel to and 



dt 



dr J &■ 



ld_ 5.^3^^ 
r dr K"^ dr 



dr 



I, (132) 



for the tilt vector. 

Comparing equations (|124|l . p25|) (present work) and (|131|) . (I132|l (Pringle, 1992), equa- 
tions H124|) and H125|l are defined for the component of angular momentum parallel to £ 
and the tilt vector, respectively. 



We may simplify the equations (|124|l . 1)1251) and (|129|) to find the general form of the 
angular momentum equation as follows 



' r Or ' 



dt 



1 d 



1 d 



-— {QiXr^nH) + -— Q^Ir^n^— + - {Q,Ir^n^£ x — ) 
1 d 



di 



dr 



1 d 



r dr 



d£^ 



dr 



1 d 



r dr 



d£ 



-- {Q\j:r'nH) + -- ( Q'^^r^n'-^ ) + ( Q'.^r^nH x - ) . (133) 



dr 



1 d 



r dr 



di" 



dr , 



